#include "misc.hpp"

#include <windows.h>
#include <GL/gl.h>
#include <GL/glu.h>
#include <GL/glut.h>
#include <stdlib.h>
#define _USE_MATH_DEFINES
#include <math.h>

#include "config.hpp"
#include "SystemMatrix.hpp"
#include "TimedVoxels.hpp"

namespace Dynamic
{
double rnd_0_1() { return (double)rand() / (double)RAND_MAX; }
double logb(double x) { return log(x) / log(logbase); }
double sqr(const double& a) { return a * a; }
double sign(double x) { if (x == 0) { return 0; } else if (x > 0) { return 1; } else { return -1; } }

void PrintTGAHeader(FILE* ofile, int width, int height)
{
	fputc(0, ofile);
	fputc(0, ofile);
	fputc(2, ofile);
	for (int i = 3; i < 12; i++) { fputc(0, ofile); }
	fputc(width % 256, ofile);
	fputc(width / 256, ofile);
	fputc(height % 256, ofile);
	fputc(height / 256, ofile);
	fputc(24, ofile);
	fputc(32, ofile);
}

void PrintTGAColor(FILE* ofile, double value, double gamma /*= 0.8*/)
{
	int r, g, b;
	if (value < 0) { value = 0; }
	if (value > 1) { value = 1; }
	//  if (value == 0) {
	//      r = g = b = 0;
	//      return;
	//  }
	/*  if (value < 0.5) {
	        r = 0;
	        g = (int)(value * 2 * 255.5);
	        b = (int)((0.5 - value) * 2* 255.5);
	    } else {
	        b = 0;
	        r = (int)((value - 0.5) * 2 * 255.5);
	        g = (int)((1 - value) * 2 * 255.5);
	    }
	*/
	r = g = b = (int)((1 - pow(value, gamma)) * 255.5);
	fputc(b, ofile); fputc(g, ofile); fputc(r, ofile);
}

void Message(char* buffer, bool up)
{
	glColor3f(1, 1, 1);
	glRasterPos2f(0.01, (up) ? 0.85 : -0.98);
	for (int i = 0; i < 80; i++) {
		if (buffer[i] == '\0') { break; }
		else { glutBitmapCharacter(GLUT_BITMAP_HELVETICA_10, buffer[i]); }
	}
}

// Gauss PDF
double Gauss2D(double dist, double sigma)
{
	if (sigma == 0) { return 0; }
	double sigma2 = sigma * sigma;
	return exp(- dist * dist / 2 / sigma2) / 2 / M_PI / sigma2;
}

// Gauss PDF
double Gauss1D(double dist, double sigma)
{
	if (sigma == 0) { return 0; }
	double sigma2 = sigma * sigma;
	return exp(- dist * dist / 2 / sigma2) / sqrt(2 * M_PI) / sigma;
}

//
// generates Poisson distributed random variables of mean lambda
//
unsigned int poisson(double lambda)
{
	double elambda = exp(-lambda);
	int n = 0;

	if (lambda < 30) {
		double p = 1;
		do {
			n++;
			p *= rnd_0_1();
		} while (p > elambda);
		n--;
	} else {
		// Rejection sampling method if lambda is too big
		double c = 0.767 - 3.36 / lambda;
		double beta = M_PI / sqrt(3.0 * lambda);
		double alpha = beta * lambda;
		double k = log(c) - lambda - log(beta);

		for (; ;)   {
			double u = rnd_0_1();
			double x = (alpha - log((1.0 - u) / u)) / beta;
			n = floor(x + 0.5);
			if (n < 0) {
				continue;
			}
			double v = rnd_0_1();
			double y = alpha - beta * x;
			double lhs = y + log(sqr(v / (1.0 + exp(y))));
			double lognfact = 0;
			for (int i = 1; i <= n; i++) { lognfact += log((double)i); }
			double rhs = k + n * log(lambda) - lognfact;
			if (lhs <= rhs) { break; }
		}
	}
	return n;
}

void PoissonRandomize(double* array, int nelems, int nsamples)
{
	double weight = 1;
	if (nsamples > 0) {
		double sum = 0;
		for (unsigned int i = 0; i < nelems; ++i) { sum += array[i]; }
		weight = sum / nsamples;
	}
	for (unsigned int i = 0; i < nelems; ++i) {
		double rvalue = (double)poisson(array[i] / weight) * weight ;
		array[i] = rvalue;
	}
}

// Base functions
double linearUniformTentBasisFunction(double dt, int timeSliceId, double t)
{
	double middle = timeSliceId * dt;
	double dist = abs(middle - t);

	if (dist > dt) {
		return 0;
	} else {
		double val =  1 - dist / dt;
		// normalize
		val /= dt; // dt is the area of the tent basis function if its height is 1.0
		return val;
	}
}

int getRandomVoxel(TimedVoxels& x)
{
	double xSum = x.Sum();
	double randVal = xSum * rnd_0_1();

	double sum = 0.0;
	int xInd = 0;
	for (; xInd < nvoxels - 1; xInd++) {
		sum += x.getAccumulatedValue(xInd);

		if (randVal < sum) {
			break;
		}
	}

	return xInd;
}

int getRandomLorForVoxel(int voxelIndex, SystemMatrix& A_LV)
{
	double alvSum = A_LV.sumByV(voxelIndex);
	double randVal = alvSum * rnd_0_1();
	double sum = 0.0;
	int l = 0;
	for (; l < nlors - 1; l++) {
		sum += A_LV.at(l, voxelIndex);

		if (randVal < sum) {
			break;
		}
	}

	return l;
}

double woodcockTrackEmissionTime(std::function<double(double, double, double)> inverseProbFunc,
								 std::function<double(double)> timeDensityFunc, double xMax, double tMax)
{
	double t0 = 0.0;
	while (true) {
		double r = rnd_0_1();
		double tT = inverseProbFunc(t0, r, xMax);
		double x_tT = timeDensityFunc(tT);

		t0 = tT;

		// Imaginary or real?
		double r2 = rnd_0_1();
		if (r2 < x_tT / xMax) {
			return tT;
		} else if (tT > tMax) {
			throw std::exception("Tmax reached");
		}

		// try again
	}
}

}